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Standard and anomalous transport in incompressible flow is investigated 
using multiscale techniques. Eddy-diffusivities emerge from the multiscale 
analysis through the solution of an auxiliary equation. From the latter it is 
derived an upper bound to eddy-diffusivities, valid for both static and time- 
dependent flow. The auxiliary problem is solved by a perturbative expansion 
in powers of the Peclet number resummed by Fade approximants and by a 
conjugate gradient method. The results are compared to numerical simulations 
of tracers dispersion for three flows having different properties of Lagrangian 
chaos. It is shown on a concrete example how the presence of anomalous 
diffusion can be revealed from the singular behaviour of the eddy-diffusivity 
at very small molecular diffusivities. 
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I. Introduction 



The problem of passive scalars diffusion in incompressible velocity fields has a theoretical 
and practical importance in many fields of science and engineering, ranging from mass 
and heat transport in geophysical flows to chemical engineering and combustion [||]. The 
main interest is in the understanding of the mechanisms leading to transport enhancement. 
Taking into account the molecular diffusion, the motion of a fluid element can be described 
by the following Langevin equation 

^ = v{x,t) + 'n{t), (1.1) 

where v{x,t) is the Eulerian incompressible velocity field at the position x and time t, 
and ?7 is a Gaussian white noise with zero mean and correlation function 

{Vi{t)vAt')) = '^Do6,,6{t-t'). (1.2) 

The coefficient Dq is the (bare) molecular diffusivity. If 6{x,t) denotes the concentration 
of tracers, the Fokker-Planck equation P| associated to (|LT1) is 

dte + {v- d)9 = Dod^e. (1.3) 



The incompressibility condition d ■ v = is explicitly used in (|1.3D . Our interest will be 



mainly concentrated on the long-time behaviour of (|I.3|) . For time scales much longer 
than the characteristic microscopic time, the evolution of 6{x,t) is dominated by long- 
wave disturbances. The equation for these slow modes can be obtained by the usual 
"hydrodynamic" analysis 

m=DS^^^{e) + ... ^,J = l,...,ci (1.4) 

where {6) is the concentration field averaged locally over a volume of linear dimensions 
much larger than the typical length I of the velocity field and d is the space dimension. 
The corrections in ( p^) involve terms containing at least three derivatives of (6), which 
can be neglected in the weak gradients limit \d{6)\/{6) <^ l^^ . Eq. ( |1.4D then reduces to 
a diffusion equation, with an effective diffusion tensor Df- (the eddy-diffusivity tensor). 
The latter has a direct practical importance since it measures the spreading for very long 
times of a spot of tracers : 

^ij = 4 ^(^»(^) " ^^»)) (^i(^) ~ ^^i))) ' ^'-^ = 1, . . . (1.5) 

■> t^oo Zt 

where x{t) is the position of a tracer at time t and the average is taken over the initial 
positions or, equivalently, over an ensemble of test particles. Note that the existence of 
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the limit in (^) ensures that the transport is a standard diffusion process, at least for 
a very large time. This is the typical situation, but there are also cases showing the 
so-called anomalous diffusion : the spreading of particles does not vary linearly with time 
but as a power law V , with 7 7^ 1 (where 7 > 1 and 7 < 1 correspond to superdiffusive 
and subdiffusive behaviours, respectively). Transport anomalies indicate the presence of 
strong correlations in the dynamics, even at large time and space scales. An interesting 
possibility is the one discussed in [||. The flow is periodic, but the Lagrangian phase 
space is a complicated self-similar structure of islands and cantori. Particles are thus 
transported in a coherent way longer and longer as Dq is deacreased, finally leading to 
anomalous diffusion in the absence of any molecular diffusion. 

The aim of this paper is using multiscale techniques to study standard and anoma- 
lous diffusion. The multiscale formalism is introduced in Section II, where the calculation 
of eddy-diffusivities is reduced to the solution of an auxiliary equation. From the latter, 
upper and lower bounds for eddy-diffusivities are derived in the general case of time- 
dependent flows. The calculation of the exact analytical expression of the eddy-diffusivity 
for parallel flows and random flows 5-correlated in time is also reviewed. Random flows 
with a short correlation time are discussed more thoroughly in Appendix 1. Numerical 
methods are generally needed to solve the auxiliary equation leading to the eddy-diffusivity 
tensor for a generic flow. Two possibilities are discussed in Section III. The first one is to 
perform a Pade resummation of the series expressing the eddy-diffusivity in powers of the 
Peclet number. The second is the use of a conjugate gradient algorithm. Both methods 
are used in Section IV to analyze three flows having a standard diffusive transport (the 
ABC 1^, H, the BC and a two dimensional time- dependent flow). The results are 
compared to numerical simulations of tracers dispersion, i.e. numerical integrations of 
([1.1|). The flows have been chosen since they can be considered as prototypes for three 
very different situations with respect to Lagrangian chaos |P, 10|, i.e. the chaotic proper- 



ties of the deterministic equation obtained from ( PTl| ) suppressing the noise. If the latter 
equation is integrable, as for the BC flow, the diffusion process is expected to be strongly 
sensitive to the detailed geometric structure of the Eulerian field and to the presence of 
molecular diffusion. For a non-integrable flow we expect a competition between coherent 
transport in the non-chaotic regions (which turns out to be dominant in the ABC flow) 
and random advection. The limiting case is the one of a strongly turbulent flow, like 
the time-dependent flow, where molecular transport can be ignored on a large range of 
scales and chaotic advection is dominant. The study of anomalous diffusion is presented 
in Section V, where the flow introduced in is analyzed. A singular behaviour of the 
eddy-diffusivity at high Peclet numbers is shown to be a signature of anomalous transport. 
A reliable procedure for predicting the presence of anomalous diffusion is thus provided. 
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II. The multiscale technique 



A general method for studying transport processes is the so-called multiscale technique 
(also known as homogenization P). The idea is to exploit the scale separation in the dy- 
namics. Specifically, let v{x, t) be an incompressible velocity field, periodic both in space 
and time. (The technique can be extended to handle the case of a random, homogeneous 
and stationary velocity field with some non-trivial modifications in the rigorous proofs of 
convergence [O].) The scalar field 6{x, t) evolves according to the Fokker-Planck equation 



([1.3|). The units are chosen in such a way that the periodicities of v are 0(1). In order 
to avoid trivial sweeping effects, the average of the velocity field over the periodicities is 
supposed to vanish. We shall be interested in the dynamics of the field 6 on large scales 
assumed to be 0(l/e), where e <C 1 is the parameter controlling the scale separation. 
Because we expect the scalar field to have a diffusive dynamics, the associated time scale 
is 0(l/e2). 

The presence of the small parameter e naturally suggests to look for a perturbative 
approach. The perturbation is however singular [|12| since a constant field is a trivial 



solution of ([TBI). The origin of this phenomenon can be grasped in the following simple 
situation. Let the large-scale field have a single wavenumber e. Because of the advection 
term in ([1.3D, a small-scale field 9 is produced and the wavenumbers spaced from those of 
V by multiples of e are generally excited. The interaction between the latter modes and 
those of V, due again to the advection term, is responsible for the transport coefficients 
renormalization. The essential shift of order e in the wavenumbers of 6 with respect 
to those of V is missed by regular perturbation expansions. Asymptotic methods, like 
multiscale techniques, are thus needed. 

In addition to the /ast variables x and t, let us then introduce s/ow variables as X = ea? 
and T = e^t. The prescription of the technique is to treat the two sets of variables as 
independent. It follows that 

d,^d, + eVr, dt^dt + e^dr, (II. 1) 

where d and V denote the derivatives with respect to fast and slow space variables, 
respectively. The solution is sought as a perturbative series 

e{x, t; X, T) = + ee(^) + e^O^^^ + . . . , (II.2) 

where the functions 9^'^^ depend a priori on both fast and slow variables. By inserting ( [II. 2|) 
and ([II. 1|) into (|]3]) and equating terms having equal powers in e, we obtain a hierarchy 
of equations. The solutions of interest to us are those having the same periodicities as 
the velocity field. The first equation, corresponding to 0{e^), is 

+ {v-d) = Do . (II.3) 
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By using Poincare inequality, one can show that for periodic solutions 

~dtj {e^'^y dv = D,J {de^'^ydv>D, {^-^^^ j (e^^y dv, (ii.4) 

where L is the spatial periodicity length of v (supposed for simplicity to be the same in 
all directions) and the integral is over the periodicity box. The inequality ( [1I.4| ) implies 
that the solution will relax to a constant with respect to fast variables, i.e. 

e^^\x,t-x,T) = e^^\x,T). (II.5) 

It can be also easily checked that the transient has no effect on the large-scale dynamics. 
The equations at order e and are 

dtO^^'^ + {v-d) e^^^ - Do d^e^^^ = -V ■ Ve(°\ (II.6) 
5^^(2) + (vd) - Do d^e^^^ = -OtO^^^ - {v ■ V)^(^) + Do V'^(°) + 2Dod ■ Ve^^\ (II.7) 
Since the equation ( [II. 6|) is linear, its solution can be written as 

e(^) {x, t;X,T)= 0^^^ (X, T) + w{x, t) ■ V0^^^ {X,T), {11.8) 

where the first term on the r.h.s. is a solution of the homogeneous equation and the vector 
field w has a vanishing average over the periodicities and satisfies 

dtw + {v- d)w-Do d'^w = -v. (11.9) 

Due to the incompressibility of the velocity field, the average over the periodicities of the 
l.h.s. in ( |ll.(j| ) and ( p.l.7| ) is zero. For the equations to have a solution, the average of 
the r.h.s. should also vanish (Fredholm alternative). The resulting solvability conditions 
provide the equations governing the large-scale dynamics, i.e. the dynamics in the slow 
variables. From (p.7|) we obtain 

c>y(^(o)) = Do V'(0(°)) - {v ■ V^W), (11.10) 

where the symbol (■) denotes the average over the periodicities. The solvability condition 
for ( [II. 6| ) is trivially satisfied, refiecting the absence of a-type effects fl^, |1^ . By plugging 
rsj) into ( [11.10[ ) we obtain the diffusion equation 

9t^(°)(X,T) = Df^ V^e^^\X,T), (11.11) 



where the eddy diffusivity tensor is 
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D§ = Do6,, - ^ [ {v,w,) + {v,Wi) ] . (11.12) 

Remark that the structure of the eddy-diffusivity tensor will reflect the rotational sym- 
metries of V and is in general non-isotropic. 
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II.. 1 Inequalities for the eddy-diffusivity 



Two important inequalities can be derived from the auxiliary equation ([II .91) and the 
expression ( |11.12| ) of the eddy-diffusivity. Let us consider the i-th and the j-th components 
of ([II .91) and multiply by wj and Wi, respectively. Taking the sum and averaging, the time 



derivative and the advective term vanish and we obtain 



^ [ {viWj) + {vjWi) ] = Do{dwi ■ dwj). (11.13) 



From (|LT3D and (|TLT|) it follows 

Df^=D,[5,^ + {dw,-dw,)]. (11.14) 

This expression of the eddy-diffusivity clearly shows that the correction to the molecular 
contribution is positive definite. Large-scale scalar transport is therefore enhanced in the 
presence of a small-scale incompressible velocity field. The cause is that for the advection- 
diffusion equation ( [r3[ ) the integral of 6*^ over the whole space is a decreasing function of 
time. When the dynamics does not possess the latter property the large-scale transport 
can actually be depleted, rather than increased. For momentum transport in Navier- 
Stokes flow the depletion can be so strong that the eddy-viscosity becomes negative, i.e. 



the average flux is in the same direction as the large-scale gradient [jT6|, |T^]. 

The second inequality also is derived from ( [11.13[ ) but it is an upper bound to eddy- 
diffusivities. Because of incompressibility, the velocity field can be expressed using a 
vector potential as v = rot A. By taking the trace of ( p.l.l3[ ) and integrating by parts, we 
obtain 

< Do{dwi ■ dwi) = -{viWi) = - {A - lotw). (IL15) 
Application of the Schwartz inequality leads to 

Do{dw,-dwi) < {A^/miotwfy/^ < {AY'^{dwi-dwiY'^, (n.l6) 

whence 

r)E //i2\ 

^< d+\-f = d+Ve\ (n.l7) 

The Peclet number is denoted by Pe. The result ( p.l.l7| ), valid for time- dependent flows 
also, generalizes a similar inequality known for time-independent velocity fields [jlT], |18[. 
The inequality ([II.17[) also provides an upper bound for each eigenvalue since the eddy- 



diffusivity tensor is positive definite. 
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II.. 2 Two exactly solvable cases 

By using multiscale techniques, the calculation of eddy diffusivities has been reduced to 
the solution of the auxiliary equation (|11.9| ). Numerical methods are generally needed to 
solve it but there are a few cases where one can obtain the solution of ( [1L9| ) analytically. 
We shall briefly review here the case of parallel flows and random flows (5-correlated in 
time. 

The peculiar property of parallel flows is that the velocity is everywhere in the same 
direction, e.g. in three dimensions 

v{x,y,z;t) = {v,{y,z;t), 0,0), (11.18) 

and Vx cannot depend on x because of incompressibility. The advective non-linearity 
V ■ dv is thus vanishing. Thanks to the latter, we can easily obtain the solution of the 
auxiliary equation ( [11. 9| ) as 

= ^^^. (11.19) 

The Fourier transforms of Vx and Wx are denoted by v and w. If F{q, u) =< \v{q, uj)\^ >, 
it follows that the eddy-diffusivity is 

Df = Djl+[ 4ttK du] ; Df = Do- (11.20) 

Here, and D\ are the components of the eddy-diffusivity tensor parallel and orthogo- 
nal to the direction of the velocity. For the Kolmogorov flow Vx = y cos y and the parallel 
eddy-diffusivity = + jW^. 

Let us now consider random flows having a short correlation time r. Neglecting the 
diffusion term in (|II.9| ) we obtain a hyperbolic equation which can be formally integrated 
along the characteristics 

w {x{a,t)]t) = — [ V {x{a, s); s) ds + w{a;0). (11-21) 
Jo 

Here, a denotes the Lagrangian initial position and the Eulerian position at time t is 

x(a,t) = a+ / V (x(a, s); s) ds. (11.22) 
Jo 

From ( [II. 21]) Taylor's expression of the eddy-diffusivity tensor immediately follows 

A?=2j( ^\M)^^'M)ds. (n.23) 
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where the Lagrangian correlation function is defined as 

rg(i-s) = {v^{x{a,t)-t)vj{x{a,sy,s)). (11.24) 

The operation (■) denotes either spatial or ensemble averaging which do coincide since the 
velocity field is supposed homogeneous, stationary and mixing. Note that the convergence 
of the integral ( [11.23| ) is not at all guaranteed. The role of a small, but non-zero, molecular 
diffusivity can be crucial in this respect In the limit where r is small, the Lagrangian 
correlation Tfj tends to the Eulerian correlation {vi{x, t) Vj{x, s)). For a signal 5-correlated 
in time 

{viix, t) vj{x, s)) = 2 Fij5it - s), (11.25) 
and the expression ( [11.23| ) reduces to 

Df< = DoS,,+F,,. (11.26) 

The corrections to this result due to a small, but finite, correlation time will be studied 
in Appendix 1. 



III. Numerical methods 



Whenever the auxiliary equation ( [11. 9|) cannot be solved exactly, numerical methods are 
needed. In this Section we shall discuss two different methods that we have used : a 
perturbative expansion and a conjugate gradient algorithm. 

In the perturbative method, the solution w of the auxiliary equation (|11.9|) is sought 
as a power series in the Peclet number Pe ~ 1/Dq: 

= Pe W + Pe^ tiJ^^) + ■ ■ • . (III.l) 

We shall concentrate on the time-independent case for simplicity. By inserting the expan- 
sion (|111.1D into ( |11.9| ) the following recursive relation is obtained 

t.« = 9-^, ^(2)^g-2^-^^^^\ ... = 9-!^:^^, . . . (III.2) 

Expressions (|III.2|) and the calculation of the average value in ( P^I.12| ) are conveniently 
handled in Fourier space, leading to 

^ = 5,, + 5:(c„),,Pe2". (III.3) 

^0 n>l 
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Here, the c„'s are numerical coefficients and the series turns out to be in Pe^, rather than 
in Pe. The contribution of order 2n + 1 in < Wj Wj > is indeed antisymmetric, as can be 
easily checked by integrating n times by parts. The series ( [III.3|) will in general converge 
for Pe < Pe* only, because of singularities in the complex plane. A reliable analytic 
continuation beyond the disc of convergence can however be performed. In |jT9[ it was 
indeed shown that the component of the eddy-diffusivity in the arbitrary direction n can 
be represented ClS db Stieltjes integral 

where ph{z) is a positive definite function, possibly singular. The poles of the eddy- 
diffusivity, considered as a function of a complex variable, are all on the imaginary axis. 
Moreover, it follows from (|III.4|) that Pade approximants of (|III.3|) have some interesting 



peculiar properties (see e.g. [0). Let us indeed denote by P^(Pe) the diagonal Pade ap- 



proximant of order n for the series ( pil.3| ) and by P^+i (Pe) the Pade approximant having 
the numerator and the denominator of degree n and n + 1, respectively. The follow- 
ing results hold for every value of the Peclet number : (i) The diagonal sequence is 
monotonically increasing and has an upper bound ; (ii) The sequence Pn+i is monotoni- 
cally decreasing and has a lower bound ; (iii) The exact value P* of the Stieltjes integral 
satisfies 

lim < < lim P" 1. (III.5) 

The difference (P^+i — P") decreases monotonically in n and provides an upper bound 
to the error due to the finite order. The quality of the resummation by a finite order 
approximant can be thus checked self-consistently. Pade approximants are very sensitive 
to the precision in the computations when the series is extended well beyond its radius 
of convergence. For small values of the molecular diffusivity, the coefficients in the series 
( [ill. 3D must be then known with very high precision. In our numerical calculations we used 



the FORTRAN multiple-precision package MP, written by R.P. Brent [^. It should be 
noted, however, that very high precision computations are quite expensive in computers 
memory costs (see next Section). 

The second method that we have used to solve the auxiliary equation (|II.9|) is a 
conjugate gradient algorithm The components of the vector w are not coupled in eq. 
([1I.9|) , which is thus equivalent to a set of scalar equations. All of them can be written in 
Fourier space as 

Aij Xj = hi z, j = 1, . . . , y. (III. 6) 

Here, V is the resolution, Xi and hi are vectors having the V components equal to the 
Fourier transform of the relevant components of w and — v, respectively. Conjugate gradi- 
ent algorithms are widely used to minimize multidimensional functions when the number 
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of dimensions is very large. The interested reader is referred to ||22| for a comparison with 



other methods (Gauss-Siedel or Minimal- Residue [PT]) in another stiff numerical prob- 



lem, the inversion of the propagator in lattice quantum chromodynamics. The solution 



of the problem ( III. 6 ) is sought by minimizing the quantity {Ax — hY over a sequence 
of directions orthogonal to the matrix A. In all the applications of the method that we 
have considered, the matrix A in (|1II.6D is sparse (quasi-diagonal). Each iteration of the 



minimization algorithm can be then performed in OiV) operations, rather than OiV"^). 
For a positive-definite matrix, the rate of convergence of the method can be shown to 
be exponential Our matrix A has actually one zero eigenvalue, corresponding to a 
constant field. The problem is nevertheless well-posed since both the velocity field and 
the solution w are orthogonal to constants, i.e. have zero average. As in any other nu- 
merical scheme, the simulations are expected to become more and more demanding as the 
molecular diffusivity becomes smaller. An increasing number of excited scales requires 
indeed a greater resolution and the rate of convergence of the method decreases when V 
and the Peclet number are increased. It is also to be checked that no eigenvalue is equal 
to zero within the numerical accuracy because of round-off errors. In the next Section 
it will turn out that the previous limitations are not very severe and do not forbid to 
perform high Peclet numbers simulations. 

We conclude this Section by briefly describing the numerical scheme used for the 
numerical simulations of tracers dispersion. The latter are done by uniformly distributing 
N particles in the periodicity box and letting them evolve according to the Langevin 
equation (pTp. The i-th diagonal element of the eddy-diffusivity tensor is then given by 

2 

(III.7) 



1 ^ 

a^it) = lim y 

> t^oo 2Nt 



1 ^ 



The indices k and j label the particles whereas the index i denotes the spatial directions 
{x,y for the two-dimensional and x,y,z for the three-dimensional case). The numerical 
integration of the Langevin equation was performed by a Runge-Kutta algorithm, modified 
to take into account the white noise term [^. The integration step was At = 0.01 and 
the total number of integration steps was 10^. This ensured a good convergence of the 
quantities ( [III. 7] ) also for the lowest molecular diffusion coefficients Dq used. The number 
of particles used was 1 000 for the three-dimensional and 2 000 for the two-dimensional 
case. 

IV. Standard diffusion 

The aim of this Section is to apply the methods previously discussed to three flows showing 
standard diffusion. The criterion in the choice of the flows is to have different mechanisms 
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of diffusion enhancement, highlighting the influence of Lagrangian chaos on transport at 
high Peclet numbers. Specifically, we have considered : 



The three dimensional ABC flow @, 0, : 

A sin(z) + C cos(?/), 

Bsm{x) + Acos{z), (IV. 1) 

C sin(y) + B cos(x). 

with A = B = C . The ABC flow is a Beltrami time-independent solution of Euler's 
equations. Eq. ( [IV. ID shows Lagrangian chaos but the phase space is also made of 
regular regions, having roughly the shape of a tube parallel to one of the three axes 
(principal vortices). 

The two dimensional BC flow 

X =Ccosiy) 
y = B cos(xj. 



obtained by projecting the flow ( |IV.1| ) onto the x — y plane and translating the 



x-coordinate by 7r/2. Eq. ( [iV.2 ) is integrable and the streamlines form a closed 



structure made of four cells in each periodicity box. 
The two dimensional time dependent flow 

X = cos{y) + sin(t/) cos(t) 



y = cos(xj + sm(x) cos(rj 

This flow is not a solution of Euler's equations anymore but it is the superposition 
of the flow ( |1V.2| ) with another flow of the same type oscillating with frequency 



uj = 1. The motivation for introducing a time dependency is to destroy all possible 
"regular islands", like the vortices in ( p.V.l| ). 



Note that both the flows ( P^V. 1|) and ( [IV. 2|) have an isotropic eddy-diffusivity tensor. 



Let us indeed consider the latter for simplicity and perform the following two operations : 
translation by vr and mirror-inversion with respect to one of the axes (e.g. x ^ n — x 
and y it — y). From the auxiliary eq. ([II. 9|) it follows that, under the previous 



operations, one of the components of w is odd and the other is even, in such a way that 
< VxWy >=< VyWx >= 0. The diagonal components are obviously equal because of the 
symmetry x ^ y. For ( [IV.l] ) the proof is similar, exploiting the fact that the group of 
symmetries of the flow is isomorphic to the cubic group |^. The flow ( [IV. 3[ ) possesses 
the symmetry x ^ y, but it is not mirror-symmetric. The diagonal components of 
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the eddy-diffusivity will then be equal but the non-diagonal component does not vanish. 
In Langevin simulations the previous symmetry properties are exploited to reduce the 
statistical fluctuations by averaging over the directions. 

In Figs. 1, 2 and 3 we present the results for the diagonal component of the eddy- 
diffusivity tensor of ( [IV. 1| ), ( [IV. 2|) and ( [IV. 3| ), respectively. The curves in each figure 
correspond to numerical simulations of the Langevin equation, the Fade method and the 
conjugate gradient algorithm. To attain the highest Feclet number the order of Fade 
approximants used is 54, 115, 29 and the number of significant digits in the computations 
is 83, 203, 40, respectively. Concerning the conjugate gradient algorithm, in Fig. 4 it 
is shown the power spectrum of the auxiliary field w for the flow ( [IV. 3[ ) at Dq = 0.01. 
It can be seen that the field is resolved enough to ensure the presence of a conspicuous 
exponentially decaying tail. The conjugate gradient algorithm turns out to be much more 
efficient at high Feclet numbers than the Fade method. The latter has the advantage of 
requiring the calculation of the coefficients of ( [111.3[ ) only : once they are computed, the 
eddy-diffusivities for all values of Dq such that the method works are available. On the 
other hand, the memory costs for high precision arithmetics are a major drawback and 
practically restrict the method to moderate Feclet numbers. 

From the high Feclet number behaviour of the eddy-diffusivities in the figures it is clear 
that the three flows have a very different dynamics. The main contribution to diffusion 
in the flow ( [IV. 1[ ) comes from the particles in the vortices, where the transport is almost 
ballistic, leading to the observed I/Dq dependence. Because of the presence of closed 
cells, a non-zero molecular diffusivity is needed to have an effective diffusion in the flow 
([IV. 2| ), as indicated by the \/D^ behaviour in Fig. 2. The transport for small molecular 
diffusivities indeed occurs by jumps from one cell to another due to the white-noise term 



in the Langevin equation ^ The probability of jumping is controlled by the 
width of boundary layers located near the separatrices and gives the square-root law. 
The flow ( [1V.3| ) is finally an example of strong Lagrangian turbulence. The particles can 
diffuse even in the absence of molecular diffusion, chaotic advection is dominant and the 
eddy-diffusivity attains a finite value indipendent of the molecular diffusivity. The figures 
show that for all the flows considered, Langevin simulations and numerical solutions of 
the auxiliary equation ([II. 9[) do agree. Moreover, in the latter method no problem of 
finite statistics and simulation times must be overcome. We conclude that multiscale 
techniques combined with an efficient numerical scheme for the solution of the auxiliary 
equation (e.g. a conjugate gradient algorithm or a pseudo-spectral code |^) provide a 
powerful tool for the calculation of eddy-diffusivities and transport properties. 
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V. Anomalous diffusion 



We shall discuss here how the multiscale formalism presented in the previous Sections 
can be used for the problem of anomalous diffusion. At a first sight it would seem that 
multiscale techniques cannot be used anymore. Consider indeed a two-dimensional static 
parallel flow ( p.l.l8| ). If the power spectrum F{q) defined in Section II. 2 is such that 

F{q) ^q"", a< 1, for g < 1, (V.l) 



then the integral in (|II.2CI| ) diverges and the eddy-diffusivity is not defined. The divergence 



is actually reflecting the fact that the transport in the direction of the flow is superdiffusive 
m, 11113,11, i.e. 

{x'{t))r^t'\ u>l/2, (V.2) 

and it is not a standard diffusion. The particles are indeed coherently swept by large- 
scale modes having wavelengths comparable (or even larger) to the typical length of the 
scalar field. Scale-separation breaks down and multiscale methods, heavily relying on this 
assumption, seem to become useless. Let us however cut out the singular part in the 
integral ( [II.20|) by defining a regularized velocity field vl, such that 




The eddy-diffusivity is now finite and exhibits a dependence 

Df{L) ~ L^-". L > 1, (V.4) 

on the cut-off length L. A standard diffusion is however observed only for spatial and 
time lengths larger than L and t* ~ L'^/Dq, respectively. For t t* the system has indeed 
p^] and for times shorter than t* it shows the same behaviour as in ( |V.2| ). By 



a crossover 



matching at t* the two different regimes, we obtain 



4 



- > 1/2. (V.5) 



For a = 0, i.e. a velocity field which is a white noise in space, ( |V.5| ) leads to u = 3/4, the 
well known result of Matheron and De Marsily ||2^ . 



In the previous example the origin of superdiffusion was related to the spatial structure 
of the velocity field. Another interesting case is the one of velocity fields with very 
long Lagrangian correlation times. The integral defining the eddy-diffusivity in Taylor's 
expression ( il.23| ) may then diverge, indicating the presence of anomalous transport for 



Dq = Q. Note however that for any Dq > Q the transport is in general a standard 
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diffusion (see (|II.20|) for an example of the role of a small molecular diffusivity) . The 
molecular diffusivity can be thus used as a regularization parameter, similarly to the 
cut-off length for parallel flows. As in the latter case, by studying the behaviour of the 
eddy-diffusivity close to the critical point (small Dq) one should be able to have some 
insights into the anomalous behaviour at the critical point [Dq = 0). Specifically, the 
examples of Section IV show that in the presence of ballistic channels the eddy-diffusivity 
varies as the inverse of Dq for small Dq, while for a system with strong Lagrangian chaos 
it tends to a constant. If denotes the eddy-diffusivity in the arbitrary direction h, we 
are thus led to interpret a small Dq behaviour 

Di ~ V, < /5 < 1 (V.6) 

as a mark of anomalous diffusion in the direction n for Dq = 0. For a practical application 
of the previous argument we have considered the velocity field [|] 

= dyip + esinz 

= —dxip + e cos z (V.7) 
= 



where 



(x + ^/2,y\ fx - V^y' 

cos X + COS + COS 



(V. 



Numerical simulations of ( |V.7|) have led the authors of [0 to conclude that the flow 
exhibits anomalous diffusion in the x — y plane for some intervals of ^-values in the range 
(0,5). In particular, 6 = 1 and e = 2.3 are such that the diffusion is standard and 
anomalous, respectively. Let us now introduce a small molecular diffusivity and calculate 
the eddy-diffusivity of the flow ( |V.7| ). The auxiliary equation is solved by the conjugate 



gradient method and the results for D^x are presented in Figs. 5 and 6. It is evident 
that for e = 1 the eddy-diffusivity tends to a constant for small Dq, while for e = 2.3 it 
is observed the behaviour D^x ~ D^^ with (3 ~ 0.7. This value is in rough agreement 
with the numerical results of using the dimensional relation /5 = 2z/ — 1 obtained 
by matching the diffusive behaviour with the anomalous behaviour (|V.2|) at the typical 



diffusive time 0(1/ Dq). The criterion ( |V.6| ) is thus confirmed and we conjecture that its 
validity is not restricted to the flow ( |V. 7[ ) only. For a generic flow, anomalies in the zero- 
diffusivity dynamics could be then captured by introducing a small molecular diffusivity 
and looking for a singular behaviour of transport coefficients. As shown in the previous 
section, the advantage with respect to simulations of the Langevin equation is that no 
problem of statistical fluctuations must be tackled. The previous procedure should then 
allow to make robust predictions on the presence of anomalous transport. 
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Appendix 1 



We shall derive here the expression of the eddy-diffusivity for incompressible flows having 
a short correlation time. Specifically, the ratio Tc/t^ between the correlation and the 
sweeping time (defined precisely later) is supposed to be small. We shall be particularly 
interested in 3D isotropic flows. By the latter we mean velocity fields invariant under 
rotations, but not in general under parity transformations. When the correlation time 
tends to zero, we recover the case of flows (5-correlated in time. For a Gaussian flow, the 
first two corrections are shown to be proportional to (tc/t^)^ : one of them is related to 
the correlation length of the flow and the second is due to helicity. The former reduces, 
while the latter increases the eddy-diffusivity, in agreement with ||33[| . 

Let v{x,t) denote a random, homogeneous and stationary incompressible velocity 
field. We shall suppose the flow to be isotropic, Gaussian and the correlation function 

{vi{x, t) v,{y, s)) = C{\t~s\) Bi,{x - y) . (A.l) 

The mean velocity is equal to zero. The temporal correlation function C{t) decays on a 
time-scale of order r^. The spatial correlation function is deflned via its Fourier transform 
as 

B,Ak) = P,Ak) ^ - -e,n h ^ , (A.2) 

where Pij = 6ij — kikj/k'^ is the solenoidal projector and eijk is the fundamental antisym- 
metric tensor. The functions E{k) and H{k) will be called the energy and the helicity 
spectrum since 

^.{v^) = C(0) / E{k) dk; {vuj) = C(0) / H{k) dk . (A.3) 



2 

The helicity is a pseudo-scalar and it is thus vanishing for flows having a center of sym- 
metry (parity-invariance). The helicity spectrum satisfles the inequality (see @]) : 

\H{k)\<2kE{k) . (A.4) 

We shall be interested in the calculation of eddy-diffusivities for very high Peclet 
numbers. The eddy-diffusivity is given in this limit by Taylor's expression ([11.231) 



2 7o 



{vi{x{a,t),t)vj{a,Qi)) + i^j 



dt . (A.5) 



Let us now suppose the correlation time Tc of the velocity fleld to be much smaller 
than the sweeping time r^. The latter is deflned as 

= (^;^ = {c{0)Jk^E{k) dk) ' ^"^'^^ 
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and it is roughly the average time it takes for a particle to travel a distance equal to 
the correlation length A. The dominant contribution in ( |A.5| ) will be given by Eulerian 
positions x{a, t) close to a : 



+ . . . 



(A.7) 



rt rt rs 

x{a,t) = a+ / dsv{a,s)+ / ds{\/iv){a,s) / ds' vi{a, s') 
Jo Jo Jo 

and the velocity v{x{a,t),t) in ( [A .51) is 

v(x(a,t),t) = v{a,t) + {Viv) {a,t) [ dsvi{a,s) + 

Jo 

{Viv){a,t) f ds {VmVi){a,s) f ds'vm{a,s') + 
Jo Jo 

-(ViV„^v){a,t) / dsvi{a,s) / ds' v^a, s') + . . . . (A.8) 
2, Jo Jo 



Eq. (|A.8|) can now be plugged into Taylor's expression ( |A.5|) , leading to 



1 



2 io 



dt 



Vi{a,t)vj{a,0)) + 
ds f ds' {{VlVmV^){a,t)vj{a,0)){vl{a,s)vm{a,s')) + 

•J 



JO 



ds I ds'{{ViVi){a,t)v.mia,s')){{VmVi)ia,s)vj{a,0)) 



(A.9) 



where homogeneity, incompressibility and the properties of Gaussian statistics have been 
exploited. Eq. ( [A.9| ) is valid for a generic Gaussian random flow and the next terms 
in the expansion are 0{tc/ts)'^. Let us now specialize (|A.9| ) to the isotropic case. The 
eddy-diffusivity tensor is then proportional to 5ij and its trace can be calculated by using 
(O) and O ■■ 



TrL) 



+ 



/roo 
E{k)dk dtC{t) - 

^ (^j E{k) dk^ (^j k^E{k) dk^ dt C{t) ds ds' C{\s-s'\) + 

i (^J H{k) dk^ dt ds C{s) ds' C{t - s') . (A.IO) 



In order to estimate the order of magnitude of the various terms in (|A.1CI|) it is convenient 
to consider the case 



1, if \t\ < 
0, otherwise. 



(A.ll) 
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When Tg — >■ a flow (5-correlated in time is obtained. In this limit, the only non- vanishing 
contribution in ( |A.10| ) is the first one, which coincides with ( [11.261 ). Both corrections in 
(|A.10|) are proportional to (tc/t^)^, as can be checked by using ( A.4|) and the Schwartz 
inequality. If the correlation function C{t) is not (|A.11|) , the constants are changed but 
not the orders of magnitude, provided the condition J C(t) dt = 1 is kept fixed. 

Note that for the correlation function (|A.11|) the helicity contribution in ( |A.10| ) is 
clearly positive while the one related to the correlation length is negative. These results 
have a simple physical interpretation. It is convenient to consider the Lagrangian corre- 
lation time which is proportional to the eddy-diffusivity. The first correction in ( [A. 101 ) 
is due to the fact that the presence of a spatial correlation length obviously reduces the 
Lagrangian correlation time. The second correction depends on the presence of helicity. 
The latter has the effect that particles move following a helix, instead of a straight line. 
The mean velocity is however the same since it depends on the energy spectrum only. 
It will then take a longer time for a particle to escape a strongly correlated region, the 
Lagrangian correlation time is longer and the eddy-diffusivity is increased. An equivalent 
remark is that a path following a helix is discriminated against tightly bending back on 
itself H. 

A Gaussian flow has been considered but the results can be easily generalized to the 
general case. If third-order moments do not vanish the first correction will in general be 
proportional to Tc/ts. We finally note that the helical term is actually the only one in 
(|A.9|) which needed to be symmetrized with respect to the indices i and j. The latter fact 
is related to Onsager's reciprocity theorem [^. When the helicity does not vanish, the 



correlation function will indeed not satisfy the time- reversibility condition Bij{x) = Bji{x) 
and the velocity field has a preferred sense of rotation The lack of parity invariance 
can be thus interpreted as a lack of symmetry with respect to time-reversal, which is 
responsible for the antisymmetry of transport coefficients. 
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FIGURE CAPTIONS 

FIGURE 1. The diagonal component Z)^ as a function of the bare molecular diffusivity 



Dq for the three dimensional ABC flow ( [IV. 1|) with A = B = C = 1. The continuous 



line is the result of the Fade method. The black and white squares are the results of the 
conjugate gradient method and direct simulations, respectively. 

FIGURE 2. The diagonal component Df^ as a function of the bare molecular diffusivity 
Dq for the two dimensional BC flow ( |1V.2| ) with B = C = 1. The continuous line is the 
result of the Fade method. The black squares are the results of the direct simulations. 
FIGURE 3. The diagonal component Df^ as a function of the bare molecular diffusivity 
Dq for the two dimensional time dependent flow ( [IV. 3| ). The continuous line is the result 
of the Fade method. The black and white squares are the results of the conjugate gradient 
method and of direct simulations, respectively. 

FIGURE 4. Log-log plot of the energy spectrum of the auxiliary field w for the ( [IV. 3[ ) 
flow at Dq = 0.01. 

FIGURE 5. The diagonal component Df^ as a function of the bare molecular diffusivity 
Dq for the three dimensional flow ( [V.7[ ), ( |V.8| ) with e = 1 (standard diffusion). 
FIGURE 6. The diagonal component Df^ as a function of the bare molecular diffusivity 
Dq for the dimensional flow ( [V.7[ ), ( [V.^ ) with e = 2.3 (anomalous diffusion). 
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